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SUMMARY 


As  part  of  Phillips  Laboratory's  Eurasian  Seismology  Program,  we  have  been  investi¬ 
gating  the  path  effects  on  seismic  waves  under  Contract  F19628-90-C-0158  during  the  period 
16  Aug  1990  *  15  Dec  1992.  Both  theoretical  and  observational  studies  were  conducted  at 
different  phases  of  this  project  to  explore  a  wide  spectrum  of  topics  that  are  directly  related  to 
the  monitoring  of  compliance  of  the  Threshold  Test  Ban  Treaty  and  the  Non  Proliferation 
Treaty: 

[1]  Quantifying  the  path  effects  on  body-wave  amplitudes  of  Novaya  Zemlya  explosions. 

[2]  Preliminary  yield  estimates  of  Novaya  Zemlya  explosions  with  path-corrected  mb . 

[3]  ISC  travel-time  inversion  in  the  Garm  Region,  Central  Asia. 

[4]  Tomographic  imaging  of  Lg  Q  structure  in  Iranian  Plateau. 

[5]  Finite-difference  modeling  study  of  ripple-fired  explosions. 

[6]  Preliminary  yield  estimate  of  Lop  Nor  explosion  of  21  May  1992  based  on  a  single¬ 
station  observation. 

The  first  two  sections  of  this  final  report  cover  the  research  performed  under  Tasks  3  and 
4,  respectively.  We  present  a  perspective  overview  of  the  whole  project  in  the  third  section, 
which  also  summarizes  the  work  done  under  Tasks  5  and  6.  The  results  obtained  under 
Tasks  1  and  2  have  been  reported  previously  in  our  first  annual  report  PL-TR-92-2042. 

The  objective  of  Task  3  is  to  obtain  the  velocity  structure  of  the  Garm  region  of  Central 
Asia.  Travel-time  residuals  of  regional  phases  recorded  at  the  IRIS  station  GAR  and  the 
CDSN  station  WMQ  were  used  to  derive  a  spherically  symmetric  structure  for  the  crust  and 
upper  mantle  in  this  region. 

Under  Task  4  we  examined  a  general  procedure  which  incorporates  the  independently 
derived  information  of  localized  path  effects  into  the  magnitude  determination.  The  goal  of  this 
exercise  is  to  quantify  the  bias  in  the  seismic  magnitude  (such  as  mb(Lg))  that  would  be 
inherent  in  a  scheme  without  fully  coupling  the  regional  propagations!  characteristics  into  the 
magnitude  determination  procedure.  Iran  was  chosen  as  a  test  case  in  this  study  because  of 
the  growing  nuclear  proliferation  concern  in  the  Middle  East.  A  tentative  zoning  partitioning 
the  Iranian  Plateau  into  six  regions  has  been  used  in  our  block  inversion  to  reveal  the  spatial 
pattern  of  the  Lg  attenuation  parameter,  Q.  Although  the  whole  Iranian  Plateau  can  be  briefly 
described  as  a  region  of  very  low  Q,  applying  a  simple  averaged  attenuation  coefficient  (Q) 
for  the  whole  plateau  would  be  inappropriate  for  mb{Lg)  calculation.  The  regional  variation  of 
the  anelastic  attenuation  parameter  is  significant  enough  that  it  must  be  taken  into  account  in 
calibrating  each  monitoring  station  for  a  reliable  magnitude  scale  for  monitoring  possible  clan¬ 
destine  tests  from  a  vast  area. 
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SECTION  I 


TRAVEL  TIME  INVERSION  IN  THE  GARM  REGION 

Christopher  $.  Lynnes'  and  Rong-Song  Jih 
Teledyne  Geotech  Alexandria  laboratory 
314  Montgomery  Street 
Alexandria,  VA  22314 


1.  INTRODUCTION 

The  propagation  of  regional  phases  is  still  problematic:  amplitudes  and  even  the  obser¬ 
vability  of  phases  can  be  highly  variable.  This  has  important  implications  for  the  application 
and  efficacy  of  earthquake/exptosion  discriminants  and  event  size  estimation.  At  least  some 
of  this  variation  is  likely  due  to  the  effects  of  lateral  heterogeneity,  particularly  in  the 
waveguides  and  near  the  boundaries  along  which  regional  phases  such  as  Pg  ,  Pn  ,  Sn  ,  and 
Lg  propagate.  The  effects  of  lateral  heterogeneity  on  regional  phase  propagation  can  be 
modeled  using  finite-difference  methods.  To  employ  these  methods,  however,  a  realistic  velo¬ 
city  model  is  required.  Block  inversion  of  travel  times  is  a  well-established  technique  for 
deriving  laterally  heterogeneous  velocity  structures.  In  order  to  simplify  the  problem  and 
improve  ray  coverage,  a  two-dimensional  geometry  is  sought:  the  high  seismicity  rate  in  the 
vicinity  of  the  IRIS  station  GAR  provides  such  an  opportunity.  The  first  step  in  this  process  is 
the  derivation  of  a  spherically  symmetric  structure  for  the  crust  and  upper  mantle  in  this 
region.  This  model  can  then  serve  as  the  starting  model  for  a  two-dimensional  velocity  inver¬ 
sion  in  the  vertical  plane. 

2.  METHODOLOGY 

The  inversion  method  used  herein  is  based  on  the  joint  inversion  of  hypocenters  and 
velocity  structure  outlined  by  Crosson  (1976).  Beginning  with  a  starting  model,  perturbations 
(Am)  to  the  previous  model  m  are  sought  to  minimize  the  travel  time  residual  vector 
d  =  jPred_Tobs,  The  model  perturbations  include  both  hypocentral  parameters  and  velocity 
structure. 

The  above  method  requires  computations  of  theoretical  travel  times,  ray  parameter  and 
9T/r)z  for  all  source-station  pairs  in  an  arbitrary  velocity  structure.  As  the  solution  of  the  two- 
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point  problem  can  be  lengthy  and  unstable  in  the  presence  of  low-velocity  zones  or  strong 
gradients,  we  use  Vidale’s  (1988)  finite-difference  technique  to  compute  the  travel  times 
throughout  a  two-dimensional  grid  for  a  source  at  the  surface.  The  various  hypocentral  partial 
derivatives  are  computed  numerically,  interpolating  from  the  grid  points  at  the  corresponding 
focal  depth  and  epicentral  distance.  Rays  are  also  traced  by  following  the  travel  time  gradient 
back  to  the  source  at  the  surface  in  order  to  compute  the  velocity  structure  derivatives. 

A  joint  hypocenter-structure  inversion  algorithm  has  been  developed  based  on  that 
described  by  Crosson  (1976).  Beginning  with  a  starting  model,  perturbations  (Arrij)  to  the  pre¬ 
vious  model  m,  are  sought  to  minimize  the  travel  t'ne  residual  vector  d(=Tjpred-Tlobs.  The  nor¬ 
mal  equations  to  be  solved  are  A(J  Amf  =  dj.  For  the  hypocentral  model  parameters,  the  com¬ 
ponents  of  A  are  the  partial  derivatives  of  the  predicted  travel  time  for  that  event,  assuming  a 
flat  earth: 

T*  -  -p  cosfo) 

Ty  -  -p  sin((j)) 

-  _ _Jl1_ _ 

vsVi-(pvT 

T,  =  1 

where  p  is  the  ray  parameter  for  a  given  source -station  pair,  b  is  the  source-station  azimuth, 
and  vs  is  the  velocity  at  the  source.  The  earth  structure  parameters,  which  are  defined  in 
terms  of  slowness  are  simply  given  by: 

TT  =  hi _ 

^Si  Vl-p2Sf2 

where  s,  is  the  slowness  (V1)  in  layer  i,  and  hj  is  the  thickness  of  layer  i  traversed  by  the  ray. 
Similar  equations  can  be  obtained  for  a  spherical  earth. 

The  resulting  matrix  A  for  inverting  p  events  and  k  layers  is: 


[H1 

0  . 

.  0 

S’l 

0 

H2 

.  0 

s2 

0 

0  . 

Hp 

SP, 

2 


where  Hj  is  the  hypocentra!  sub-matrix  and  S,  is  the  velocity  structure  sub-matrix  for  event  i: 


T*.1 

Ty,i 

ta« 

TU 

T*,2 

V 

Tz,2 

Tx,N, 

Ty,Ni 

Tz,N, 

Tt,N, 

Ts,t 

H 

v> 

••  Tski 

T, 

T 

T 

S12 

S2S 

S|(2 

Ts,N, 

TSjn, 

"  TskN, 

Nj  is  the  number  of  observations  for  event  i. 

In  reality,  each  row  of  A  and  d  is  weighted  according  to  the  precision  and  quality  (impul¬ 
sive  or  emergent)  of  the  arrival  time  report  contributing  to  d. 

The  equation  AAm  =  d  is  then  solved  for  Am  by  multiplying  both  sides  by  AT: 

AtA  Am  =  ATd 


(AtA)-'  AtA  Am  =  (AtA)_1  ATd 

The  matrix  ATA  is  a  sparse  matrix  ot  doubly-bordered  block-diagonal  form.  The  inversion 
thus  lends  itself  well  to  conjugate  gradient  techniques,  such  as  that  described  by  Press  et  al. 
(1986).  Unfortunately,  the  matrix  ATA  is  often  close  to  being  singular,  particularly  since  there 
are  strong  tradeoffs  between  the  velocity  in  the  source  layer  and  the  partial  derivative  with 
respect  to  depth.  The  near-singularity  renders  the  matrix  inversion  unstable.  To  get  around 
this  problem,  we  decided  to  adopt  the  Sherman-Morrison  inversion  procedure  (Press  et  al. . 
1986),  in  which  the  matrix  inverse  is  found  for  the  purely  block-diagonal  matrix,  and  then  suc¬ 
cessively  modified  for  eaJi  row  and  column  making  up  the  borders.  Increasing  the  number  of 
events  thus  results  in  an  increase  of  only  order  N2  for  the  inversion  of  the  normal  equations. 

In  the  actual  joint  inversion  algorithm,  the  model  perturbations  are  computed  from  inver¬ 
sion  using  a  starting  model,  and  are  then  added  to  the  starting  model.  The  procedure  contin¬ 
ues  until  model  perturbations  are  sufficiently  small  for  a  given  iteration,  or  the  improvement  of 
the  travel  time  residual  vector  is  sufficiently  small.  At  each  iteration,  both  the  travel  time  and 
ray  parameter  must  be  recomputed  for  all  source-station  pairs  for  the  current  structure.  In 
order  to  solve  the  two-point  problem  for  a  large  number  of  pairs,  we  use  Vidale's  (1988) 
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algorithm  for  computing  travel  times  everywhere  within  a  grid  using  finite-difference.  At  each 
iteration,  this  grid  need  be  computed  only  once  for  a  source  at  the  surface;  travel  times  and 
ray  parameters  are  then  obtained  from  the  grid  point  at  the  appropriate  hypocentral  depth  and 
epicentral  distance. 

The  implementation  of  the  algorithm  comprises  three  programs: 
bvnormeqn : 

constructs  the  normal  equations. 
hv  invert'. 

inverts  for  the  solution,  i.e. ,  the  model  perturbations. 
hvadjust: 

adds  the  model  perturbations  to  the  previous  model. 

This  separation  of  the  inversion  program  from  the  others  allows  us  to  vary  parameters 
with  each  iteration,  such  as  the  damping  added  to  the  matrix  ATA,  as  well  as  to  recover  easily 
from  inversion  failures  due  to  non-convergence  of  the  conjugate  gradient  method.  The  joint 
inversion  algorithm  has  been  tested  on  synthetic  cases  ano  is  now  being  applied  to  the  ISC 
arrivals. 

3.  DATA  AND  RESULTS 

Hypocenters  and  arrivals  have  been  obtained  from  International  Seismologicai  Center 
(ISC)  tapes  for  the  years  1971-1982.  They  have  been  dearchived  and  reformatted  for  use  by 
software  for  earth-structure  inversion.  The  travel  time  data  actually  used  in  this  study  consist 
of  those  recorded  at  stations  less  than  1000  km  (9°)  away.  The  ISC  locations  (Figure  1)  are 
used  as  the  initial  guesses  for  hypocentral  parameters.  Events  were  selected  carefully  to 
include  close  stations  so  that  the  tradeoff  between  depth,  origin  time  and  source  layer  velocity 
could  be  resolved. 

Several  studies  {e.g. ,  Carter  ef  a/.,  1991)  indicate  the  crust  to  be  very  thick  in  this 
region.  Figure  2  shows  reduced  travel  time  residuals  for  events  with  ISC  depths  o*  10-15  km. 
The  travel  time  curve  for  the  Soviet  GAR  model  (which  has  a  mantle  velocity  of  about  8.22 
km/s)  is  superimposed.  The  match  is  fairly  good  at  this  depth,  with  the  thick  crust  confirmed 
by  the  large  crustal-mantle  phase  crossover  distance.  Figure  3  shows  the  residuals  for 
events  with  ISC  depths  between  40  and  45  km.  The  deterioration  at  depth  may  indicate  that 
ISC  depths  are  systematically  overestimated  due  to  the  velocity  model. 
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A  starting  model  similar  to  the  GAR  model  was  used  in  the  inversion.  The  mantle  velo 
city  was  fixed  at  8  2  km/s  in  order  to  stabilize  the  layer  directly  above.  However,  this  velocity 
appears  to  be  well  constrained  by  the  travel  time  slooe  from  400-1000  km  (Figures  2  3).  The 
crustal  depth  was  fixed  at  65  km.  Figure  4  shows  the  results  of  the  joint  inversion.  The 
resulting  model  is  similar  to  the  GAR  model.  The  velocity  gradient  is  slightly  steeper  in  the 
upper  crust.  Below  this  is  a  35-km  thick  lower  crust  of  about  7  km  s,  slightly  lower  than  the 
GAR  mode!  velocity.  The  relocated  depths  are  plotted  against  the  ISC  depths  in  Figure  5.  As 
indicated  in  Figure  3,  the  relocated  depths  are  generally  less  than  the  ISC  depths. 

4.  CONCLUSIONS 

A  joint  inversion  algorithm  for  hypocenters  ana  velocity  structure  has  been  developed 
based  on  the  algorithm  of  Crosson  (1976).  Vidale's  (1988)  finite  difference  travel  time  algo¬ 
rithm  and  the  Sherman-Morrison  inversion  of  sparse  matrices  have  been  used  to  speed  exe¬ 
cution  for  large  numbers  of  events.  A  joint  inversion  for  hypocenters  and  velocity  structure  in 
Central  Asia  yields  a  1 -dimensional  model  with  a  thick  crust  consisting  a  30-km  crust  with  a 
steep  gradient  and  a  roughly  homogeneous  35-km  tower  crust  of  about  7  km/s.  Focal  depths 
appear  to  be  systematically  overestimated  by  the  ISC.  This  model  can  be  used  as  the  start 
ir.g  model  for  a  2-dimensional  tomographic  inversion  of  a  profile  in  central  Asia. 
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Figure  2.  Reduced  travel  times  tor  events  10-15  km  in  the  shaded  region  of  Figure  1  with  ISC  depths 
of  10-15  km.  (The  inversion  events  in  this  depth  range  are  a  subset  of  this  set).  The  predicted  travel 
time  curve  for  the  Soviet  GAR  mode!  reported  by  Carter  el  al.  (1991)  is  superimposed. 
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Figure  3.  Same  as  Figure  2  for  a  depth  range  of  40-45  km  The  small  reduced  travel  times  for  close 
stations  suggests  the  ISC  depths  are  systematically  overestimated 
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1.  INTRODUCTION 

An  accurate  characterization  ot  the  propagation  of  regional  phases  due  to  explosion 
sources  in  regions  of  high  proliferation  concern  (such  as  Iran)  requires  fairly  accurate  crustal 
models.  Iran  lies  on  the  Alpine-Himalayan  seismo-tectonic  belt,  within  a  wide  band  of  seismic 
activity  connecting  the  Hindu  Kush  and  eastern  Mediterranean  highly  seismic  regions.  The 
Zagros  Main  Thrust  Fault  constitutes  the  boundary  ;etween  the  Iranian  and  Arabian  Plates. 
Central  Iran,  northwest  Iran,  and  the  south  Caspian  areas  have  low  levels  of  seismicity  com¬ 
pared  to  those  of  the  narrow  belts  surrounding  them  (McEvilly  and  Razani,  1971;  Nowroozi, 
1971,  1972,  1976;  McKenzie,  1972;  Dewey  and  Grantz,  1973;  Hedayati  et  al.,  1976;  Chandra 
et  al.,  1979;  Quinttmeyer  and  Jacob,  1979;  Rowshandel  et  al.,  1981;  Ambrasseys  and  Mel¬ 
ville,  1982;  Jackson  and  McKenzie,  1984,  1988;  and  many  others).  This  is  regarded  as  a 
highly  heterogeneous  and  complicated  geological  environment  (Leith,  1992). 

The  seismic  structure  of  the  Iranian  Plateau  has  been  discussed  in  several  studies.  For 
instance,  using  travel  time  residuals  derived  from  data  recorded  at  WWSSN  stations  in  Iran 
(MSH,  SHI,  TAB),  Chen  et  al.  (1980)  inferred  an  uppermost  mantle  P-wave  velocity  of 
8.010.1  km/s  for  Iran.  They  also  suggested  that  the  crust  dips  toward  the  south-southeast  at 
about  1°.  Thus  if  the  crustal  thickness  were  34  km  in  the  north,  it  would  reach  about  49  km  in 
the  south.  Kadinsky-Cade  et  al.  (1981)  also  derived  the  average  upper-mantle  velocity 
beneath  Iran  with  the  same  single-station  method.  Asudeh  (1982)  applied  the  two-station 
method  to  earthquakes  along  the  interstation  line,  and  he  derived  a  P- wave  velocity  of 
7.7010.05  km/s  with  a  crust  43  km  thick  between  southern  Iran  and  the  SRO  station  MAIO. 
This  profile  crosses  the  Tabas  area.  Asudeh  (1982)  also  derived  simple  velocity  profiles  for 
several  other  regions  with  MAIO  and  ILPA  (Iranian  Long  Period  Array)  data.  Asudeh  (1982) 
suggests  that  the  crust  to  the  north,  east,  and  central  parts  of  the  Iranian  Plateau  is  approxi¬ 
mately  10  km  thicker  than  at  its  western  extreme.  Roughly,  the  crust  thickens  from  SSW  to 
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NNE  across  the  Zagros  Main  Thrust  Fault,  contrary  to  the  assertion  made  by  Chen  et  al. 
(1980).  Wallace's  (1991)  over-simplified  two-block  model  of  Iran  seems  to  support  Asudeh's 
(1982)  interpretation.  A  point  that  should  be  made  here  is  that,  although  a  few  crustal  models 
have  been  proposed  during  the  past  three  decades,  the  requirement  to  identify  small  seismic 
events  has  raised  a  need  to  refine  the  available  crustal  models  for  more  detailed  and  more 
accurate  regionalized  profiles  based  on  newly  available  data  sources  and  methodologies. 

This  study  covers  two  topics.  In  the  first  part,  we  establish  a  regionalized  Lg  Q  map  for 
Iranian  Plateau  with  a  7-zone  tomographic  inversion.  In  the  second  part,  we  illustrate  the 
importance  of  applying  the  2-dimensional  Lg  Q  map  to  mb{Lg)  calculation.  The  bias  in 
mb{Lg)  due  to  erroneous  constant  Q  model  is  quantified.  Due  to  the  preliminary  nature  of 
the  result,  however,  it  is  not  our  intention  to  explain  the  cause  of  the  regional  variation  of  Lg 
Q  in  Iranian  Plateau.  Mitchell  and  Hwang  (1987)  find  that  the  variation  of  Lg  Q  values 
between  Applachian  and  Rocky  Mountains  can  easily  be  produced  by  accumulations  of  sand 
stone  and  shale  of  Mesozoic  age  and  younger.  By  contrast,  throughout  most  of  the  Western 
United  States,  neither  the  low  Lg  Q  values  which  have  been  observed  nor  the  regional  varia¬ 
tion  of  those  values  can  be  explained  by  accumulation  of  low -Q  sediments.  Mitchell  and 
Hwang  (1987)  suggest  that  the  explanation  requires  low  and  laterally  varying  values  of  O  in 
the  crystalline  crust  for  Western  United  States.  To  explore  the  cause  of  the  regional  variation 
of  Lg  Q  in  Iranian  Plateau  would  require  more  data  and  further  investigation. 
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2.  NUTTU'S  Lg  ATTENUATION  MEASUREMENTS 

Nuttli  (1980)  selected  99  Iranian  earthquakes  which  occurred  during  1972-1974.  For 
these  earthquakes,  the  body-wave  magnitude.  mb ,  varied  from  3.7  to  6.0,  with  a  median  of 
4.8,  and  the  focal  depth  varied  from  0  to  94  km,  with  a  median  of  43  km.  Nuttli  measured  the 
7  of  109  paths  recorded  at  three  WWSSN  [World  Wide  Standard  Seismograph  Network]  sta¬ 
tions  in  Iran:  Mashad  (MSH),  Tabriz  (TAB),  and  Shiraz  (SHI).  However,  Nuttli  (1980)  did  not 
tabulate  the  events  he  used.  Instead  of  cross  referencing  the  epicentral  coordinates  published 
in  the  Bulletin  of  the  International  Seismological  Centre,  we  digitized  Nuttli's  hand-drawn  maps 
to  retrieve  his  epicenters.  Figure  1  shows  the  re-constructed  109  Lg  paths. 

Nuttli  (1980)  measured  the  ", sustained  maximum  motion",  namely  the  amplitude  equaled 
or  exceeded  by  the  three  largest  amplitude  waves,  of  the  vertical-component  Lg  waves  with 
period  around  1  second.  A  "normalized”  amplitude.  A*,  is  defined  as  follows:  the  amplitude 
readings,  A(a)  (in  pm),  are  first  corrected  for  the  effects  of  geometrical  spreading  and  disper¬ 
sion  with  the  formula  appropriate  for  an  Airy  phase  ( cf .  Ewing  et  at.,  1957);  then  they  are 
normalized  according  to  Nuttli’s  empirical  scaling  that  an  Iranian  event  with  ISC  mb  5.0  would 
have  a  hypothetical  Lg  amplitude  of  270  pm  at  10-km  distance.  All  events  are  thus  “equal¬ 
ized"  to  magnitude  5.0: 

A* (A)  s  A(A)  •  A0  333  •  [sin(A/111.1)}°5  •  10—^  ,  [1] 

If  there  is  no  anelastic  attenuation  or  other  scattering  mechanism,  the  normalized  amplitude, 
A*(A),  should  be  1  at  any  distance,  A.  Otherwise,  A* (A)  would  be  smaller  than  1,  and  it  is 
defined  as  exp[-y(A-1 0)],  where  y  is  the  anelastic  absorption  coefficient.  When  the  normal¬ 
ized  amplitudes  A’ (A)  are  plotted  as  a  function  of  A  on  a  semi-logarithmic  scale,  the  slope  of 
the  straight  line  connecting  [10km,  0]  and  [A,  ln(A*(A))]  would  be  exactly  -y.  Nuttli  used  this 
time-domain  approach  to  determine  the  attenuation  coefficient,  y,  of  109  paths  in  his  study 
(cf.  Figure  3  of  Nuttli,  1980).  The  average  y  of  these  109  paths  is  0.00442  km-1,  very  close 
to  that  of  0.0048  for  coastal  California  derived  by  Herrmann  (1980).  Solid  and  dashed  lines  in 
Figure  1  represent  those  paths  with  y  smaller  and  larger,  respectively,  than  the  average. 
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43.000  64.000 


Figure  1.  109  Lg  paths  recorded  at  3  WWSSN  stations  in  Iran  (MSH,  TAB,  and  SHI)  during  1972-1976 
for  which  the  attenuation  parameter,  y,  was  readily  measured  by  Nuttli  (1980)  Solid  and  dashed  lines 
represent  those  with  y  value  smaller  and  larger  than  the  average  (0.0044  km  '),  respectively. 
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3.  TOMOGRAPHIC  IMAGING  METHOD 


The  datum  for  seismic  tomographic  imaging  is  the  integral  of  some  physical  parameter 
along  a  specified  path  through  the  medium.  For  example,  the  travel  time  accumulated  along 
a  ray  path  between  a  source  and  a  receiver  can  be  expressed  as  the  integrated  slowness;  the 
amplitude  decay  can  be  expressed  as  the  integrated  attenuation.  The  application  of  tomo¬ 
graphic  imaging  to  construct  the  2-dimensional  map  of  lateral  heterogeneity  in  the  crust  is 
essentiaiiy  the  same  as  the  method  that  has  been  widely  used  for  decades  in  X-ray  tomogra¬ 
phy  or  radioastronomy. 

Suppose  there  are  M  seismic  sources  recorded  at  some  or  all  of  N  stations,  and  sup 
pose  that  these  events  and  stations  are  spread  over  K  seismotectonic  regions.  For  coding 
and  computational  simplicity,  the  region  boundaries  were  chosen  on  a  0.333°  by  0.333°  grid. 
Thus  each  "seismotectonic  province"  is  composed  of  many  square  cells  of  size  0.333°  by 
0.333°.  We  assume  that  the  attenuation,  yk,  varies  from  one  province  (region)  to  another.  The 
total  path  attenuation  can  be  decomposed  as  the  sum  of  the  attenuation  incurred  in  each 
region  that  the  ray  has  traversed.  That  is, 

K 

ZTk  Ak(i,j)  =  7(i.j)  A(i.j)  .  for  i=1 . M  ;  j=1 . N  [2] 

k=l 

where  y(i.j)  is  the  path  attenuation  that  Nuttli  (1980)  already  measured,  and  A(i,j)  is  that  epi- 
central  distance  from  the  i’th  event  to  the  j'th  station.  On  the  left-hand  side,  yk  is  the  (unk¬ 
nown)  attenuation  coefficient  in  the  k  th  region  and  Ak(i,j)  is  the  distance  traversed  in  the  k’th 
region  along  the  path  from  the  i'th  event  to  the  j’th  station.  Ak(i,j)  is  computed  as  the  sum  of  all 
segments  of  the  great-circle  path  of  the  (i.j)'th  ray  that  fail  in  a  square  cell  comprising  the  k’th 
region.  Thus  Ak(i,j)  would  be  0  if  the  (i.j)’th  ray  does  not  cross  any  square  cel!  of  the  k'th 

region.  The  determination  of  Ak(i,  j)  for  i=1 M,  j=1 N,  is  the  crucial  step  of  this  linear 

tomographic  imaging  problem.  Once  the  elements  of  matrix  [A^i,  j)]MxN  are  known,  it  is 

straightforward  to  solve  for  the  unknown  vector,  (y,,  y2 . yK]',  with  the  standard  least-squares 

techniques. 
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4.  REGIONAL  VARIATION  OF  Q  IN  IRANIAN  PLATEAU 


There  have  been  different  opinions  regarding  how  Iranian  Piateau  should  be  partitioned 
in  terms  of  the  seismotectonic  properties.  Nowroozi  (1976)  presented  23  seismotectonic  pro¬ 
vinces  for  Iran  based  on  638  relocated  earthquakes  as  well  as  24  instrumentally  located  epi¬ 
centers  given  by  others.  Berberian  (1979)  largely  disagreed  with  the  boundaries  Nowroozi 
(1976)  proposed,  and  he  suggested  that  the  only  justifiable  partitioning  of  the  Iranian  Plateau 
should  be:  [1]  two  Marginal  Fold  Belts  of  Zagros  (in  southwestern  Iran)  and  Kopeh  Dagh  (in 
northeastern  Iran)  resting  on  the  Arabian  and  Turan  platforms,  respectively;  [2]  the  Central 
Iranian  ranges  and  basins  between  them;  and  finally,  [3]  the  post-colored  melange-orhiolitic 
flysch  belt  of  east  and  southeast  Iran.  Berberian’s  (1979)  characterization  of  Iranian  Plateau 
implies  that  some  seismotectonic  provinces  in  Nowroozi's  (1976)  map  need  to  be  merged. 

If  the  ray  paths  sample  each  portion  of  the  region  of  interest  fairly  well,  then  the  tomo¬ 
graphic  imaging  can  be  utilized  to  reveal  the  regional  characteristics  without  the  need  to 
impose  any  specific  a  priori  zoning.  Thus  the  tomographic  inversion  technique  can  be  used 
to  test  various  hypotheses  regarding  the  partitioning  of  Iran,  whenever  the  data  resolution  per¬ 
mits.  If,  however,  the  data  coverage  is  poor  in  a  certain  portion,  as  in  the  case  of  Nuttli's 
(1980)  data  set,  then  applying  some  regionalization  based  on  the  seismotectonic  information 
is  necessary  and  reasonable.  Figure  2  shows  a  tentative  regionalization  which  partitions  Iran 
into  6  regions:  Zagros  Range,  Lut  Block,  East  Iran  Range,  Central  Iran  Range, 
Elburz/Caspian  region,  and  Great  Kavir/Esfaphan/Rezaiyeh. 

The  tomographic  imaging  result  with  this  partitioning  is  listed  in  Table  1.  East  Iran 
Range  and  western  Afghanistan  have  the  poorest  ray  coverage  in  Nuttli's  data  set,  and  hence 
the  uncertainty  in  the  associated  Q  estimate  is  high.  For  other  regions,  the  standard  error  in 
the  Q  value  is  around  10%.  Both  the  Zagrcs  Range  and  the  Lut  Block  show  a  large  y  of 
0.005  km-1,  roughly  corresponding  to  a  Q  of  1 81  ±1 2  and  183±18,  respectively.  The  Kopet 
Dagh,  Shahrud  Doruneh,  and  the  Qom  region  also  seem  to  have  a  Q  slightly  smaller  than  the 
average.  On  the  other  hand,  the  Elburz  Province  and  central  Iran  have  a  O  of  about  250  for 
1  Hz  Lg  waves,  which  may  be  the  highest  O  value  in  Iran.  Despite  the  large  uncertainty  in 
the  O  value  for  the  western  Afghanistan  area,  the  estimated  O  of  201  is  in  agreement  with 
Nuttli's  (1981)  earlier  work  in  which  he  pointed  out  that  Afghanistan,  northern  Pakistan,  and 
northern  India  have  crustal  Q  values  comparable  to  those  of  Iran  and  the  western  United 
States. 
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Figure  2.  A  preliminary  regionalization  which  partitions  Iran  into  6  seismotectonic  regions  Zagros 
Range,  Lut  Block,  Eastern  Iran  Range,  Central  Iran  Range,  Elburz/Caspian  region,  and  Great 
Kavir/EsfapharVRezaiyeh.  Block  inversion  indicates  that  both  the  Zagros  Range  and  the  Lut  Block  show 
a  large  y  of  0.005  km-1,  roughly  corresponding  to  a  O  ol  181  +  12  and  183+18.  respectively  The  Kopet 
Dagh,  Shahrud  Doruneh,  and  the  Corn  region  also  seem  to  have  a  O  slightly  smaller  than  the  average 
The  Elburz  Province  and  central  Iran  have  a  O  ol  about  250  tor  1  Hz  Lg  waves,  which  may  be  the 
highest  O  value  in  Iran. 
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Table  1.  Regional  yand  Q  of  1  Hz  Lg  in  Iran 


Region 

y  ('/km) 

Oo 

Central  Iran  Plateau 

0.003710.0006 

246138 

Elburz  +  Caspian 

0.003510.0004 

257126 

East  Iran  Range 

0.003810.0013 

248191 

Western  Afghanistan 

0.004510.0025 

2011112 

Great  Kavir  +  Esfaphan  +  Rezaiyeh 

0.004310.0003 

2091.6 

Zagros  Range 

0.004910.0003 

181112 

Lut  Block 

0.GO4 910.0005 

183118 

Table  2.  Crustal  O  Results  of  Various  Regions* 


Author 

Geometrical 
Spreading  + 
Dispersion 

f  (Hz) 

Region 

Method  and 
Distance  Range 

Nuttli  (1973) 

1500 

a-S/6 

1 

ENA 

time  domain  amplitudes 

Street  (1976) 

900 

a-5/6 

1 

US  &  Canada 

E  of  Rockies 

lime  domain  amplitudes 

A  -  200-4500  km 

Horner  et  al.  (1978) 

1500 

A  ~5'6 

1 

Canadian  shield 

time  domain  amplitudes 

A  -  300-7000  km 

Herrmann  (1980) 

229 

1 

Berkeley,  CA 

j  Mechler  et  al.  (1980) 

[ 

1 

400  f  0  44 

A-5/6 

0.5-8 

France 

A  -  350-1100  km 

time  domain  max  amps 
bandpass  filtered  data. 

Mitchell  (1981) 

-900  f  0  2 

0.25-1 

ENA 

A  -  500-2000  km 

theoretical  modeling 
of  spectral  shape 

Nuttli  (1981) 

800  f  0  5 

1500  f 0  2 

-1  -5 

ENA 

CUS 

A  <  1000  km 

Nuttli  (1981) 

200  f  0  7 

-1-5 

S.  Calif. 

Lg  coda,  A  <  1000  km 

1  Dwyer  et  al.  (1983) 

1 

1280  f  04 

a-5/6 

1-10 

Central  US 

time  domain  max  amps 
from  bandpass  filtered 
data,  A  =  200-2000  km 

'l  partially  adaDted  Pom  Gone*  et  al.  (1987). 
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Table  2.  Crustal  Q  Results  of  Various  Regions*  (Continued) 


Author 

Q0  $ 

Geometrical 
Spreading  + 
Dispersion 

f  (Hz) 

Region 

Method  and 
Distance  Range 

Singh  &  Herrmann  (1983) 

900  f  0  35 

1000  f  01 

1200  f  0  2 

900  f  0  25 

A'56 

0. 5-3.5 

NENA 

S.  Appalachians 
Central  US 

Central  Central  US 

Lg  coda,  300-800  km 

Pulti  (1984) 

660  f  0  4 

0.75-10 

NENA 

coda  amps,  -  50-600  km 

Chen  et  al.  (1984) 

(400-500)  f  0  0 
(120-200)  f0  2 

-1 

-1 

Beijing  area 
Yun-Nan,  China 

coda  amps, 
coda  amps. 

Raoof  &  Nuttli  (1985) 

(130-350)  f<0  4'0  7) 
(420-580)  f<0  2_0  3) 
(580  -980)  ft°-°-°-2> 

A'5'6 

0.4-1. 4 

Andes 

Patagonian 

ES/ 

Hasegawa  (1985) 

900  f  0  2 

A-172 

0.6-20 

Canadian  shield 

spectra  of  ground  accei  , 
A  -  70-900  km 

Shin  (1985) 

(500±50)  f  (0-6_0  7) 
(500±50)  f  (°  6~°  7) 

A- 5/6 
a-1/2 

1-10 

1-10 

Eastern  Canada 
Eastern  Canada 

time  domain  maw  amps 
freq  domain  amp  spect-a 

A  -  100-1000  km 

Campillo  et  al.  (1985) 

290  f  0  52 

a-5/6 

0.5-10 

France 

A  =  150-2000  km 

time  domain  n  ax  amps 
&  theoretical  modeling  j 

Chavez  8  Priestley  (1986) 

206  f  0  68 

214(±50)  f0  54±0  09 

A"1'2 

0.3-10 

0.3-5 

US  Great  Basin 
explosions 

US  Great  Basin 
earthquakes 

frequency  domain  amps, 

A  «  200-500  km 
"  ; 

*<  j 

Gupta  &  McLaughlin  (1987) 

800  f  0  32 

400  at  1  Hz 

1100  at  1  Hz 

1400  at  1  Hz 

A-5/6 

0.5-7 

1 

EUS 

CUS 

EUS-CUS 

Appalachians 

amplitudes  from  PSRV 

A  <  1100  km 

Gonez  et  at.  (1987) 

1000  f  0  35 

A-s/e 

0.5-13 

Central  &  ENA 

time  domain  max  amps 

A  -  500-3000  km  | 

Chun  et  al.  (1987) 

1100  f 0  19 

A-1'2 

0.6-10 

Eastern  Canada 

spectral  ratios  i 

A  »  53-210  km 

Xie  8  Mitchell  (1990) 

low 

A~5/e 

E.  Africa  Rift 

, 

Chun  et  al.  (1992) 

564 (±53)  f°  49±0  06 
283(±5)  f°  72±0  02 

a-5/6 

0.1-1.0 
0.3-3. 5 

Asia 

Beijing  area 

CDSN  +  IRIS  I 

A  <  1000  km 

•)  Partially  adaptod  from  Gonez  ot  ai  (1987). 
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5.  IMPLICATION  FOR  SEISMIC  MONITORING 


Attenuation  of  short-period  Lg  -wave  motions  is  of  concern  for  a  number  of  seismological 
problems.  For  engineering  purposes,  Lg  is  often  the  wave  group  responsible  for  the  largest 
and  most  prolonged  ground  motion  observed  on  accelerograms,  and  hence  it  determines  the 
"design  motion"  at  a  site.  For  seismicity  studies,  Lg  amplitudes  can  be  used  to  estimate  1- 
Hz  mb  values  for  small  earthquakes  for  which  the  P  wave  cannot  be  seen  at  teleseismic  dis¬ 
tances  (Nuttli,  1973).  For  discrimination  between  earthquakes  and  explosions  at  regional  dis¬ 
tances,  the  ratio  of  Lg  to  P  amplitudes  has  been  investigated  in  many  studies  (e.g. ,  Bland- 
ford,  1981;  Lynnes  et  at.,  1990;  Lynnes  and  Baumstark,  1991). 

For  the  purpose  of  calculating  mb(Lg)  in  Iran,  Nuttli’s  original  Lg  formula  can  be  rewrit¬ 
ten  in  an  equivalent  form: 

mtiLg)  .  3.6372  +  logA(A)  ♦  IfcgW  ♦  >1  +  ■  [3] 

where  y  =  —  j-  -  ,  Q(f)  =  O0  '  t  . 

A  is  the  epicentral  distance  in  km,  A(A)  is  the  observed  Lg  amplitude  measured  in  the  time 
domain  in  urn  [microns]  at  the  epicentral  distance  of  A  km.  For  instance,  a  seismic  source  in 
Iran  with  1-sec  Lg  amplitude  of  270  urn  at  10  km  epicentral  distance  would  correspond  to  a 
mb{Lg)  of  3.6372  +  2.4314  +  0.3333  -  1.4019  +  0.0000  =  5.000,  same  as  what  Nuttli’s  (1986) 
formula  would  give.  Nuttli's  (1986)  later  study  of  Lg  suggested  that  a  different  constant  of 
4.0272  is  required  for  Eastern  North  America  and  Central  Asia  regions  (Jih,  1992).  The  bias 
of  0.39  m.u.  could  reflect  the  difference  in  Lg  excitation  relative  to  mb  (ISC)  between  Eastern 
U.S.  (or  Central  Asia)  and  Iran.  Alternatively,  it  could  reflect  the  mb  bias  between  these  two 
areas.  This  issue  will  be  discussed  later. 

Once  the  regionalized  y  map  is  available,  an  individual  path  y  for  an  arbitrary  source- 
station  pair  can  simply  be  computed  as  the  weighted  sum  of  the  y^'s  of  the  subregions  that 
the  ray  path  traverses: 

K  Ak 

y=ZYk-  r-.  [4] 

k=1  A 

It  is  clear  from  [3]  that  an  erroneous  path  y  would  yield  a  mb(Lg)  bias  which  increases  with 
the  distance.  Furthermore,  this  error  is  independent  of  the  actual  source  size  or  the  quality  of 
the  amplitude  measured  at  the  recording  station.  It  is  the  bias  solely  due  to  inaccurate  calibra¬ 
tion  of  the  propagation  effect. 
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Figure  3  shows  the  spatial  pattern  of  apparent  mb(Lg )  corrections  that  would  be  needed 
at  WWSSN  station  SHI  and  IRIS  station  ASH  (Ashkhabad),  respectively,  if  a  constant  y  of 
0.0044  km-1  (viz. ,  the  average  y  across  Iran)  was  assumed  in  computing  mb(Lg)  at  these  two 
stations.  For  each  hypothetical  epicenter,  the  mb(Lg)  correction  required  is  independent  of 
the  amplitude  actually  observed  ( cf .  Equation  [3]). 

Partitioning  Iran  into  21  regions  as  proposed  by  Nowroozi  (1976)  (Figure  4)  in  our  block 
inversion  would  yield  more  prominent  spatial  variations  in  the  resulting  Q  map  (Figure  5),  and 
hence  a  stronger  variation  in  mb(Lg)  residual  as  well.  The  regional  variation  of  the  anelastic 
attenuation  parameter  in  Figure  5  is  significant  enough  that  it  needs  to  be  taken  into  account 
in  calibrating  each  monitoring  station  for  a  reliable  magnitude  scale  in  monitoring  possible 
clandestine  tests  from  a  vast  area.  Unless  this  has  been  done,  adding  more  stations/arrays 
to  the  monitoring  network  may  not  provide  substantial  improvement  in  reducing  the  error  in 
the  network-averaged  magnitude  based  on  regional  phases.  In  fact,  since  the  bias  increases 
with  distance,  the  bias  in  mb(Lg)  would  be  much  smaller  if  the  magnitude  were  based  on  the 
nearest  station  alone  (as  compared  to  using  a  network  of  sparse  stations  for  simple  averag¬ 
ing).  Figure  6  shows  the  spatial  pattern  of  apparent  mb(Lg)  corrections  as  a  function  of 
hypothetical  hypocenters  that  would  be  needed  for  a  network  of  3  stations  (SHI,  TAB,  and 
MSH)  if  a  constant  y  of  0.0044  km-1  were  assumed  in  computing  mb(Lg)  at  these  three  sta¬ 
tions.  The  residuals  are  determined  with  respect  to  7-  (top)  and  21 -zone  (bottom)  tomo¬ 
graphic  inversions,  respectively.  Note  that  magnitudes  based  on  the  network  average  is  not 
necessarily  better  than  those  based  on  the  nearest  station  alone,  if  the  path  attenuation  is  not 
carefully  accounted  for.  This  should  be  obvious  by  comparing  Figures  3  and  6. 
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Figure  3.  The  spatial  pattern  ot  apparent  mb(l,])  corrections  as  a  (unction  ot  hypothetical  hypocenter:. 
tnat  would  be  needed  at  WWSSN  station  SHI  (top)  and  IRIS  station  ASH  (bottom)  it  a  constant  y  ot 
0  0044  km'’  {viz,  the  average  y  across  Iran)  were  assumed  in  computing  mhU  0)  at  these  two  stations 
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Figure  4.  A  more  detailed  partitioning  of  Iran  as  proposed  by  Nowroozi  (1976)  which  partitions  Iran  into 
21  regions 
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Figure  5.  Partitioning  Iran  into  21  regions  as  proposed  by  Nowroozi  (1976)  in  our  block  inversion  would 
yield  more  prominent  spatial  variations  in  the  resulting  O  map  as  well  as  a  more  complicated  pattern  of 
mb(L;j)  residual. 
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Figure  6.  The  spatial  pattern  ol  apparent  mh{T5}  corrections  as  a  (unction  of  hypothetical  hypocenters 
that  would  be  needed  lor  a  network  ol  3  stations  (SHI,  TAB,  and  MSH)  if  a  constant  y  ol  0  0044  km'1 
were  assumed  in  computing  mb{Lg )  at  these  three  stations  Top  and  bottom  show  the  results  of  7-  and 
21 -zone  tomographic  inversion,  respectively  Magnitude  based  on  the  network  average  is  not  neces¬ 
sarily  better  than  that  based  on  the  nearest  station  alone,  il  the  path  attenuation  is  not  carefully 
accounted  for 


27 


6.  CONCLUSION 


Although  the  whole  Iranian  Plateau  can  be  briefly  described  as  a  region  of  very  low  Q, 
applying  a  simple  averaged  attenuation  coefficient  (Q)  for  the  whole  plateau  would  be  inap¬ 
propriate  for  measuring  mb(Lg).  The  regional  variation  of  the  anelastic  attenuation  parameter 
is  significant  enough  that  it  needs  to  be  taken  into  account  in  calibrating  each  monitoring  sta¬ 
tion  for  a  reliable  magnitude  scale  in  moni^ring  possible  clandestine  tests  from  a  vast  area. 
Unless  this  has  been  done,  adding  more  ations/arrays  to  the  monitoring  network  may  not 
provide  substantial  improvement  in  reducing  the  error  in  the  network-averaged  magnitude 
based  on  regional  phases.  In  fact,  since  the  bias  increases  with  distance,  the  bias  in  mb(Lg) 
would  be  much  smaller  if  the  magnitude  were  based  on  the  nearest  station  alone  (as  com¬ 
pared  to  using  a  network  of  sparse  stations  for  simple  averaging).  This  situation  may  be  very 
different  from  that  of  *  .  •  coring  a  specific  nuclear  test  site  for  which  the  empirical  site- 
dependent  corrects  ,  ‘  .n  can  be  applied  afterwards  to  the  station  magnitudes  even  when  the 
path  Os  are  not  fully  known  beforehand.  For  istance,  Jih  (1992)  recently  demonstrated  a  joint 
inversion  pm^edure  which  can  determine  the  Lg  -based  explosion  size  estimate  as  well  as  the 
path  attenuation  coefficient  simultaneously.  Such  a  technique  would  be  most  useful  for  cali¬ 
brating  the  regions  where  underground  nuclear  tests  or  other  man-made  seismic  events  fre¬ 
quently  take  place.  For  regions  like  Iran  where  relatively  very  little  is  known  about  what  the 
explosion  signature  should  look  iike,  however,  a  regionalized  y  map  (e.g. ,  Singh  and 
Herrmann,  1983,  for  North  America;  and  that  in  Kadinsky-Cade  et  a/.,  1982,  for  Sn  )  for  each 
crustal  phase  of  interest  should  be  established  and  applied  to  the  routine  magnitude  computa¬ 
tion  procedure. 

7.  RECOMMENDATIONS 

Nuttli  (1980)  suggested  that  a  seismic  source  in  Iran  with  the  ISC  bulletin  mb  5.0  should 
excite  Lg  amplitude  of  270  microns  at  a  10-km  extrapolated  distance,  whereas  the  seismic 
sources  in  Central  Asia  and  eastern  North  America  with  comparable  ISC  bulletin  mb  would 
excite  a  Lg  amplitude  of  110-150  microns  (Nuttli,  1986).  Possible  explanations  for  this 
apparent  bias  oi  0.26-0.39  m.u.  include  [1]  uncertainty  in  the  ISC  bulletin  mb  values  (Nuttli, 
1981),  [2]  differences  in  Lg  excitation  relative  to  mb  ,  or  [3j  differences  in  the  upper  mantle 
property  which  could  cause  a  bias  in  mb  .  Many  Lg  studies  have  used  the  ISC  mb  for  "nor¬ 
malizing"  their  Lg  magnitude  scale  {e.g.,  page  2146  of  Nuttli,  1986a;  page  128  of  Israelson, 
1992).  Without  careful  re-processing,  the  ISC  bulletin  mb  values  typically  would  be  associ¬ 
ated  with  large  uncertainty.  Therefore,  it  is  very  important  to  adopt  a  better  mb  in  calibrating 
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the  crustal  phases  in  Iran  (or  any  other  area  of  proliferation  concern)  in  the  first  place.  Magni¬ 
tudes  based  on  the  short-period  teleseismic  recordings  corrected  for  the  station  amplifications 
as  well  as  the  path  effects  (due  to  the  focusing/defocusing  near  the  source  region  etc.)  such 
as  those  in  Jih  and  Wagner  (1992ab)  may  better  serve  such  a  purpose. 

Wallace  (1991)  reports  that  the  central  Iran  region  has  more  efficient  Lg  propagation 
than  the  surrounding  area,  qualitatively  in  agreement  with  our  result  in  Table  1.  Many  previ¬ 
ous  studies  {e.g.,  Ruzaikan  eta/.,  1977;  Kadinsky-Cade  eta!.,  1981;  Ni  and  Barazangi, 
1982;  Pomeroy  et  at..  1982;  Wallace,  1991)  have  pointed  out  that  Lg  ,  Sn  ,  and  P„  do  not 
necessarily  have  the  same  propagation  characteristics.  For  instance,  the  inefficient  Sn  propa¬ 
gation  beneath  the  northern  Iranian  Plateau  identified  by  Kadinsky-Cade  et  at.  (1981)  seems 
to  permit  efficient  Lg  propagation  (Figure  6).  The  spectral  and  phase  discrimination  metho¬ 
dologies  should  take  this  into  account.  A  regionalized  y  map  {e.g.,  Figure  5  for  Lg  in  Iran; 
the  map  derived  by  Singh  and  Herrmann,  1983,  tor  North  America;  and  that  in  Kadinsky-Cade 
et  at.,  1982,  for  Sn )  for  each  crustal  phase  of  interest  should  be  established  and  applied  to 
the  routine  magnitude  computation  procedure  such  as  the  one  in  the  knowledge-based 
automatic  processing  system  (Bratt  et  at. ,  1991). 

Fault-plane  solutions  for  many  suitably  located  earthquakes  in  the  Middle  East  have 
been  determined  by  McKenzie  (1972)  and  Jackson  and  McKenzie  (1984).  We  recommend 
assembling  events  recorded  at  regional  as  well  as  at  teleseismic  distances  according  to  the 
clustering  of  earthquakes  with  similar  focal  mechanism.  Preferably  for  each  crustal  phase  at 
least  10  events  in  each  clustered  area  ( cf .  Figure  7)  should  be  measured.  The  redundancy  of 
data  in  each  group  will  provide  better  constraints  for  inferring  the  empirical  station  and 
region-specific  path  corrections  (Jih  and  Wagner,  I992ab;  Jih.  1992),  and  it  will  thereby  result 
in  more  precise  magnitude  scales  for  both  regional  and  teleseismic  phases.  The  determina¬ 
tion  of  seismic  magnitudes  for  events  occurring  in  isolated  suspected  spots  will  also  benefit 
from  this  exercise  significantly,  since  more  accurate  station  corrections  will  be  available. 
Although  several  crustal  models  of  the  Middle  East  have  been  suggested  during  the  past 
three  decades,  the  task  of  identifying  small  seismic  events  requires  a  refinement  of  the  avail¬ 
able  models  for  more  detailed  and  more  accurate  regionalized  profiles.  Consequently,  more 
tomographic  inversions  of  the  Middle  East  would  be  very  useful. 

To  be  more  specific,  we  recommend  first  to  assemble  broadband  digital  data  of  recent 
seismic  events  of  Iran  (and  Iraq)  recorded  at  IRIS  stations/arrays  (ASH,  KIV,  and  GNI)  and 
then  to  complement  this  data  set  with  digital  data  of  MAIO  (SRO)  and  analog  recordings  of 
historical  events  recorded  at  WWSSN  and  other  former  Soviet  stations  nearby.  We  suggest 


that  phase  parameters  of  Lg ,  Rayleigh  waves,  and  other  regional  phases  be  measured,  and 
then  inverted  for  the  regionalized  anelastic  attenuation  coefficients  as  well  as  for  the  crustal 
structures  of  Iran  {and  Iraq).  Our  preliminary  result  along  these  lines  as  presented  in  this 
study  can  be  improved  further  by  using  an  enlarged  data  set  along  with  a  more  carefully 
defined  seismotectonic  map.  We  recommend  that  the  mb  (of  relatively  larger  events)  and 
mb(Lg)  be  determined  with  two  simultaneous  inversion  procedures  which  have  been  tested 
with  Soviet  events  (Jih  and  Wagner,  1992ab;  Jih,  1992).  The  resulting  refined  mb  values 
should  then  be  applied  to  normalize  the  regional  and  local  magnitudes  of  Iranian  events.  More 
detailed  phase  discrimination  studies,  such  as  those  in  Lynnes  and  Baumstark  (1991)  and 
Lynnes  ef  al.  (1990),  should  be  conducted  for  the  Iranian  region  with  the  improved  phase 
parameters  to  be  obtained  in  the  manner  we  suggest. 
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STUDIES  OF  SEISMIC  WAVE  PROPAGATION  IN  EURASIA 

Rong-Song  Jih  Christopher  S.  Lynnes 

Teledyne  Geolech  Hughes  STX 

314  Montgomery  Street  7601  Ora  Glen  Drive,  Suite  300 

Alexandria,  VA  22314  Greenbelt,  MD  20770 

1.  CONTRACT  NO-:  F19628-90-C-0158  (16  Aug  1990  -  15  Dec  1992) 

2.  MOST  IMPORTANT  SCIENTIFIC  RESULTS 

Under  this  Air  Force-funded  contract,  we  have  been  conducting  both  theoretical  and 
observational  studies  covering  a  wide  spectrum  of  topics  directly  related  to  monitoring  the 
compliance  of  both  the  Threshold  Test  Ban  Treaty  [TTBT]  and  the  Non-Proliferation  Treaty 
[NPT].  Each  of  the  following  paragraphs  summarizes  a  paper  that  has  been  or  will  be 
presented/published.  The  investigator  who  performed  the  specific  task  is  also  identified. 

2.1  ISC  Travel  Time  Inversion  in  the  Garm  Region,  Central  Asia1  (Lynnes) 

The  propagation  of  regional  phases  is  still  problematic:  amplitudes  and  even  the  obser¬ 
vability  of  phases  can  be  highly  variable.  This  has  important  implications  for  the  application 
and  efficacy  of  earthquake/explosion  discriminants  and  event  size  estimation.  At  least  some 
of  this  variation  is  likely  due  to  the  effects  of  lateral  heterogeneity,  particularly  in  the 
waveguides  and  near  the  boundaries  along  which  regional  phases  such  as  Pg  ,  Pn  .  Sn  ,  and 
Lg  propagate.  The  effects  of  lateral  heterogeneity  on  regional  phase  propagation  can  be 
modeled  using  finite-difference  methods. 

In  order  to  employ  these  methods,  a  realistic  velocity  model  is  required.  Block  inversion 
of  travel  times  is  a  well-established  technique  for  deriving  laterally  heterogeneous  velocity 
structures.  In  order  to  simplify  the  problem  and  improve  ray  coverage,  a  two-dimensional 
geometry  is  sought:  the  high  seismicity  rate  along  the  great  circle  arc  between  the  IRIS  sta¬ 
tion  GAR  and  and  the  CDSN  station  WMQ  provides  such  an  opportunity.  Travel  times  for 
regional  phases  are  used  to  derive  a  spherically  symmetric  structure  for  the  crust  and  upper 
mantle  in  this  region.  This  model  is  then  starting  model  for  a  two-dimensional  velocity  inver¬ 
sion  in  the  vertical  plane. 

'Reported  in  Semi-annual  RSD  Status  Reports  #1  and  «2  daled  14  February  and  19  August.  1991,  respectively  See  also 
Section  I  ot  this  report 
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2.2  Path  Effects  on  Body-wave  Amplitudes  of  Novaya  Zemlya  Explosions2  (Jih) 

The  standard  procedure  used  in  estimating  the  source  size  of  underground  nuclear 
explosions  using  mb  measurements  has  been  to  separate  the  station  terms  from  the 
network-averaged  source  terms.  The  station  terms  thus  derived  actually  reflect  the  combina¬ 
tion  of  the  path  effect  and  the  station  effect,  when  only  those  events  in  a  close  proximity  are 
utilized.  If  worldwide  explosions  are  used  in  the  inversion,  then  the  path  effect  tends  to  be 
averaged  out  at  each  station.  In  either  case,  the  effect  due  to  the  propagation  path  alone 
would  not  be  obvious. 

The  major  research  goal  of  this  contract  is  to  improve  our  fundamental  understanding  of 
the  seismic  wave  propagation  in  Eurasia.  To  achieve  this  goal,  we  further  decompose  the 
station  amplification  effect  with  a  joint  inversion  scheme  that  simultaneously  determines  the 
seismic  source  size,  the  path  terms,  and  the  receiver  terms.  Short-period  P- wave  amplitudes 
of  217  worldwide  underground  nuclear  explosions,  including  28  blasts  from  Novaya  Zemlya, 
recorded  at  118  WWSSN  stations  have  been  used  in  one  single  inversion  to  isolate  the  pro¬ 
pagation  complexities  affecting  the  P- wave  amplitudes.  For  all  28  Novaya  Zemlya  events  in 
our  WWSSN  database,  the  new  mb  factoring  procedure  provides  more  stable  mb  measure¬ 
ments  across  the  whole  recording  network  with  a  reduction  in  the  fluctuational  variation  by  a 
factor  of  up  to  3.  A  typical  reduction  factor  in  the  variation  is  about  2.0,  and  the  the  worst  case 
is  1.4. 

The  inferred  path  terms  are  then  compared  against  the  travel-time  residuals  to  character¬ 
ize  the  propagation  paths.  Our  result  indicates  that  paths  from  the  northern  test  site  in 
Novaya  Zemlya  to  stations  in  North  America  have  systematically  faster  arrivals  and  smaller 
amplitudes,  suggesting  a  profound  defocusing  effect  on  the  first  arrivals;  while  stations  in  Ire¬ 
land,  Scotland,  Spain,  Bangladesh,  northern  India,  Pakistan,  Korea,  and  Kenya  report  slow 
arrivals  and  large  amplitudes,  suggesting  a  focusing  effect.  Amplitudes  for  paths  to  Green¬ 
land,  Iceland,  Alaska,  Turkey,  Germany,  Luzon,  Zimbabwe,  Italy,  Pueto  Rico,  Ethiopia,  and 
Hawaii,  however,  seem  to  be  controlled  by  the  anelastic  attenuation  with  slow  rays  also  asso¬ 
ciated  with  small  amplitudes,  and  fast  rays  associated  with  large  amplitudes.  A  strong  corre¬ 
lation  between  P- wave  amplitude  and  Lg  detection  at  teleseismic  distance  is  also  observed. 


?Reportod  in  Semi  annual  R&D  Status  Report  #3  dated  16  February.  1992,  and  Scientific  Report  #1.  PL-TR  92-2042 
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2.3  Yield  Estimates  of  Novaya  Zemlya  Explosions:  Preliminary  Result3  (Jih) 

As  a  byproduct  of  the  study  described  in  Section  2.2,  we  have  derived  the  yield  esti¬ 
mates  of  Novaya  Zemlya  explosions  based  on  the  path-corrected  mb  values.  Assuming  the 
basic  coupling  and  the  mantle  condition  at  Novaya  Zemlya  are  comparable  to  those  at 
Eastern  Kazakhstan,  the  mb  bias  relative  to  NTS  at  50  KT  level  using  the  path-corrected 
mb(Pmax )  values  is  inferred  as  0.25  and  0.36  magnitude  unit  for  Novaya  Zemlya  and  Semi- 
palatinsk,  respectively.  The  mb(Pmax)  bias  of  0.11  between  Semipalatinsk  and  Novaya  Zem¬ 
lya  could  be  largely  due  to  the  difference  in  pP  interference  between  these  two  test  sites 
(rather  than  the  seismic  coupling).  The  relative  source  size  determined  by  Burger  et  at. 
(1986)  and  the  theoretical  log('F^)  :yield  scaling  are  combined  to  extrapolate  our  mb  scaling 
to  the  higher  end.  The  resulting  yield  estimates  of  Novaya  Zemlya  explosions  range  from  2  to 
2100  KT,  with  peak  values  at  550  KT  and  65  KT  for  events  before  and  after  1976,  respec¬ 
tively,  which  are  in  reasonable  agreement  with  those  in  previous  studies. 

Also  included  in  the  Annual  Scientific  Report  No.  1  is  a  complete  listing  of  path-corrected 
mb  values  and  yield  estimates  of  Semipalatinsk  explosions  which  are  used  as  our  baseline  in 
calibrating  Novaya  Zemlya  explosions.  First  motion  of  the  initial  short-period  P  waves 
appears  to  be  a  very  favorable  source  measure  for  explosions  fired  in  hard  rock  sites  under¬ 
lain  by  the  stable  mantle  (such  as  Semipalatinsk). 

2.4  Regional  Lg  Q  Variation  In  the  Iranian  Plateau  Revisited4  (Jih) 

We  examined  a  general  procedure  which  incorporates  the  independently-derived  infor¬ 
mation  of  localized  path  effects  into  the  magnitude  determination.  The  goal  of  this  exercise  is 
to  quantify  the  bias  in  the  seismic  magnitude  (such  as  mb[Lg))  that  would  be  inherent  in  a 
scheme  without  fully  coupling  the  regional  propagational  characteristics  into  the  magnitude 
determination  procedure.  Iran  was  chosen  as  a  test  case  in  this  study  because  of  the  growing 
nuclear  proliferation  concern  in  the  Middle  East. 

A  tentative  zoning  partitioning  Iranian  Plateau  into  six  regions  {viz. ,  the  Zagros  Range, 
the  Lut  Block,  East  Iran  Range,  Central  Iran  Range,  the  Elburz/Caspian  area,  and  the  Great 
Kavir/Esfaphan/Rezaiyeh  region)  has  been  used  in  our  block  inversion  to  reveal  the  spatial 
pattern  of  the  Lg  attenuation  parameter,  O.  We  used  109  observations  recorded  at  3 
WWSSN  stations  in  Iran  (MSH,  TAB,  and  SHI)  for  which  the  y  was  readily  measured  by  Nuttli 

■’Reported  in  Somi  annual  R&O  Status  Roport  »3  dated  16  Fobruary,  1992,  and  Scientific  Report  #1.  PL-TR-92  2042 
“Reported  tn  Semi  annual  R&D  Status  Roport  #4  dated  14  August,  1992  See  also  Section  II  of  this  report 
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(1980).  The  weighted  average  y  of  these  paths  is  0.0044  km  \  very  close  to  the  0.0048  for 
coastal  California  derived  by  Herrmann  (1980).  It  is  also  in  agreement  with  the  Lg  coda  Q 
map  recently  compiled  by  Xie  and  Mitchell  (1992)  for  Eurasia.  In  our  inversion  with  a  more 
detailed  regionalization,  both  the  Zagros  Range  and  the  Lut  Block  show  a  large  y  of  0.005 
km-1,  roughly  corresponding  to  a  Q  of  181±12  and  183±18.  respectively.  The  Kopet  Dagh, 
Shahrud  Doruneh,  and  the  Qom  region  also  seem  to  have  a  Q  slightly  smaller  than  the  aver¬ 
age.  On  the  other  hand,  the  Elburz  Province  and  central  Iran  have  a  O  of  about  250  for  1  Hz 
Lg  waves,  which  could  be  the  highest  O  value  in  Iran. 

Although  the  whole  Iranian  Plateau  can  be  briefly  described  as  a  region  of  very  low  Q, 
applying  a  simple  averaged  attenuation  coefficient  ( Q )  for  the  whole  plateau  would  be  inap¬ 
propriate.  The  regional  variation  of  the  anelastic  attenuation  parameter  is  significant  enough 
that  it  must  be  taken  into  account  in  calibrating  each  monitoring  station  for  a  reliable  magni¬ 
tude  scale  in  monitoring  possible  clandestine  tests  from  a  vast  area.  Unless  this  has  been 
done,  adding  more  stations/arrays  to  the  monitoring  network  may  not  provide  substantial 
improvement  in  reducing  the  error  in  the  network- averaged  magnitude  based  on  regional 
phases.  In  fact,  the  bias  in  mb(Lg)  would  be  much  smaller  if  the  magnitude  was  solely  based 
on  the  nearest  station  alone  (as  compared  to  using  a  network  of  sparse  stations  for  simple 
averaging).  This  situation  is  very  different  from  that  of  monitoring  a  specific  nuclear  test  site 
for  which  the  empirical  site-dependent  correction  term  can  be  applied  afterwards  to  the  station 
magnitudes  even  when  the  path  Os  are  not  fully  known  in  advance. 

2.5  Theoretical  Investigation  of  Ripple-fired  Explosions5  (Jih) 

A  major  issue  for  the  Non-Proliferation  Treaty  is  the  discrimination  of  large  chemical 
explosions  from  possible  clandestine  or  small  nuclear  tests.  Unless  discrimination  is  possible, 
the  numerous  mining  blasts  could  give  ample  opportunity  for  concealing  clandestine  tests. 
Rippled-fired  explosions  are  commonly  used  for  fragmenting  rocks  during  quarry  and  open-pit 
mining.  The  periodicity  inherent  in  the  ripple  firing  could  produce  a  seismic  reinforcement  at 
the  frequency  of  the  delay  between  shots  or  rows.  Since  ripple  firing  is  essentially  a 
sequence  of  single  explosions  arranged  and  fired  using  a  series  of  time-delayed  detonations 
between  adjacent  shot  holes  or  rows  of  shot  holes,  it  has  been  suggested  that  the  convolution 
of  a  single  explosion  with  a  comb  function  of  variable  spacing  and  variable  amplitude  can  be 
used  to  study  the  distinctive  signature  of  ripple  firing  (Stump,  1988;  Smith,  1989).  Baumgardt 
ana  Ziegler  (1988)  delicately  demonstrated  that  the  incoherent  array-stack  spectra  can  be 

''Reported  in  Semi  annual  R&D  Status  Report  #4  dated  14  August  1992 
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used  to  identify  some  multiple  shots  recorded  at  NORSAR.  By  superpositioning  the  waveform 
due  to  a  single  shot  with  proper  time  delay,  they  were  able  to  model  the  source  multiplicity 
under  the  assumption  that  the  spatial  spreading  of  the  shots  are  negligible  with  respect  to  the 
traveling  distance. 

There  are,  however,  wave  excitation  characteristics  of  ripple-fired  explosions  which  are 
not  predicted  by  such  spectral  or  waveform  superposition  approaches.  Figure  1  of  our  4th 
semi-annual  R&D  status  report  shows  the  vertical  component  of  the  displacement  wave  field 
of  7  ripple-fired  explosions  detonated  with  a  rupture  velocity  of  2.5  km/s  in  a  granitic  half 
space  Neither  Rf)  nor  the  S‘  was  emitted  in  the  single  explosion  case  using  the  same  shot 
depth.  However,  a  very  clear  Rg  is  excited  in  the  forward  direction  for  the  ripple  firing.  In 
fact,  one  of  the  standard  industrial  practices  in  reducing  potential  damage  caused  by  ground 
vibration  is  to  detonate  ripple  shots  in  a  direction  away  from  the  buildings  (Dick  et  at.,  1983). 
An  immediate  implication  of  this  exercise  on  the  discrimination  problem  is  that  the  lack  of  Rg 
is  not  necessarily  indicative  of  deep  sources.  Furthermore,  although  path  effects  such  as 
anelastic  attenuation  and  scattering  by  shallow  heterogeneity  and  topography  in  the  upper 
crust  can  reduce  Rg  significantly  as  demonstrated  by  Jih  et  al.  (1988)  and  McLaughlin  and 
Jih  (1986,  1987),  the  reported  lack  of  Rg  in  many  seismograms  from  known  quarry  blasts 
could  also  be  due  to  the  intrinsic  source  effect  (such  as  the  scooting  pattern)  rather  than  the 
path  effects. 

Several  previous  theoretical  studies  with  spectral  superposition  approach  (e.g. ,  Smith, 
1989;  Greenhalgh,  1980)  suggest  that  for  ripple  firings  with  simple  configurations,  the  spec¬ 
trum  will  be  reinforced  at  every  1/At  Hz,  where  the  At  is  the  delay  time  between  the  shots. 
However,  our  numerical  experiments  indicate  that  even  with  a  very  simple  linearly  distributed 
shot  array  in  the  half  space,  the  frequencies  at  which  the  amplitude  reinforcement  could  occur 
also  depend  on  the  relative  azimuth  angle  and  take-off  angle  with  respect  to  the  rupture  direc¬ 
tion.  Figure  2  of  our  4th  semi-annual  R&D  status  report  gives  the  spectral  ratios  of  P- wave 
synthetics  to  that  due  to  a  single  explosion.  It  is  clear  that  reinforcement  does  occur  at  exactly 
every  50  Hz  along  the  orthogonal  direction,  consistent  with  Smith's  prediction  (1987,  1989). 
For  other  directions,  however,  the  frequency  at  which  the  enhancement  would  occur  is  clearly 
off  due  to  the  Doppler  Effect.  This  phenomenon  suggests  that  the  shift  of  the  spectral  rein¬ 
forcement,  whenever  observable,  could  be  used  as  a  simple  discriminant  between  multiple 
explosions  and  a  single  shot.  It  could  also  provide  an  alternative  explanation  of  why  the  spec¬ 
tral  reinforcement  was  not  observed  at  the  expected  frequency  at  RSON  for  the  mining  blast 
C  detonated  at  Mesabi  Range,  August  6  1986  (Smith,  1989).  Kim  et  al.  (1991)  report  that 
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seismograms  from  the  same  quarry  recorded  at  the  same  station  could  have  distinct  spectral 
characteristics.  They  attribute  this  phenomenon  to  the  considerable  fluctuation  in  delay  times 
and  sub-charge  sizes.  The  numerical  modeling  study  we  conducted  provides  another  possible 
explanation,  however. 

2.6  Lop  Nor  Explosion  of  21  May  1992  as  Observed  at  ANP  (Anpu,  WWSSN),  Taiwan 

(Jih) 

The  "megaton"-level  explosion  detonated  on  21  May  1992  at  the  test  site  between  Lop 
Nor  and  Bosten  Lake,  Xinjiang,  is  believed  to  be  the  largest  Chinese  underground  nuclear 
test  conducted  so  far.  The  raw  body-wave  magnitudes  (viz. .  log(ATT)  +  B(A))  measured  at 
WWSSN  station  ANP  (Anpu,  121.517°E,  25.183°N,  A=3523km)  are  5.72.  6.17,  and  6.43, 
based  on  the  first  motion  (zero-to-peak  or  Pa  ),  the  “b"  phase  (first  peak-to-first  trough  or  Pb  ), 
and  the  largest  peak-to-peak  (Pmax)  amplitudes,  respectively.  Applying  the  empirical  station 
correction  of  0.143  magnitude  unit  [m.u.]  along  with  the  path  correction  of  -0.045  m.u.  derived 
by  Jih  and  Wagner  (1992a)  to  bring  the  station  recordings  at  ANP  in  alignment  with  the 
optimal  worldwide  network  average,  the  final  station  mb  values  at  ANP  are  set  at  5.82  (Pa  ), 
6.27  (Pb  ),  and  6.52  (Pmax),  respectively.  Since  there  is  no  evidence  of  a  low  velocity  zone 
[LVZ]  in  the  area  of  eastern  Tien  Shan  Mountains  or  the  test  site  (Matzko,  1992),  it  would 
seem  reasonable  that  the  upper-mantle  conditions  at  Lop  Nor  test  site  are  very  similar  to 
those  at  Eastern  Kazakh  where  the  major  nuclear  test  sites  of  the  former  Soviet  Union  are 
located.  As  a  result,  the  magnitude-yield  calibration  curve  derived  from  the  Pa  phases  of 
Eastern  Kazakh  explosions  ( cf .  Jih  and  Wagner,  1992a,b)  could  be  adequate  for  Lop  Nor 
explosions  as  well  The  calibration  curves  for  Pb  and  Pmnx  phases  can  then  be  inferred  from 
ihe  empirical  Pb  -Pa  and  Pmax -Pa  regressions,  respectively.  Under  the  assumption  that  the 
source  coupling  as  well  as  the  anelastic  attenuation  at  Lop  Nor  test  site  are  comparable  to 
those  at  Eastern  Kazakh,  and  further  that  the  corner-frequency  effect  is  negligible,  the  yield  of 
21  May  1992  explosion  is  estimated  as  350  KT  based  on  the  mb  values  measured  at  ANP 
alone.  If  we  assume  that  Lop  Nor  area  has  a  t*  intermediate  between  those  at  Nevada  test 
site  (U.S.A.)  and  Eastern  Kazakh,  then  our  yield  estimate  would  become  550  KT  instead.  This 
seems,  however,  very  unlikely  as  the  alternative  assumption  lacks  geophysical  support.  The 
news  media  report  that  the  nuclear  device  was  emplaced  in  a  shaft  somewhat  over  900 
meters  deep  which,  if  correct,  would  suggest  a  yield  most  likely  in  the  range  of  400-480  KT 
and  no  more  than  550  KT.  Thus  350  KT  and  550  KT  can  be  regarded  as  the  lower  and 
upper  bounds,  respectively,  of  our  preliminary  yield  estimate. 
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